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We study rectification of a current of particles moving in a spatially pe- 
riodic potential under the influence of time-periodic forces with zero mean 
value. If certain time-space symmetries are broken a non-zero directed cur- 
rent of particles is possible. We investigate this phenomenon in the framework 
of the kinetic Boltzmann equation. We find that the attractor of the Boltz- 
mann equation completely reflects the symmetries of the original one-particle 
equation of motion. Especially, we analyse the limits of weak and strong 
relaxation. The dc current increases by several orders of magnitude with 

decreasing dissipation. 
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In recent years, a lot of efforts have been devoted to the study of a directed current of 
particles in a spatially periodic potential under the simultaneous action of an external time- 
dependent field with zero mean []J. Most of the previously done work was connected with 
transport in the presence of a stochastic external field (noise) . Another set of problems 
concerns systems driven by a deterministic periodic force |5|-[7] with applications to current 
generation in semiconductor superlattices ||, a free particle moving in a non- Newtonian 
liquid ||, two-dimensional conducting electron gas in arrays of triangular shaped quantum 
dots jltj, to name a few. For more information we refer the interested reader to |TTJ. 

The model Hamiltonian of a particle of unit mass and its dynamical equation of motion 
can be written as: 



H = — + U(x)-xE(t) , x = -U\x)+E(t) . (1) 

Here U is a spatially periodic potential U(x) = U(x + 2tt) , E is an ac field E(t) = E(t + T) 
with zero mean value and frequency uj = 2ir/T. 

A recent approach to the problem of current rectification [Q was based on the analysis of 
symmetries of the dynamical equations of motion with and without additional dissipation. 
Consider trajectories generated from an ensemble of initial conditions in phase space. A 
dc current can be calculated by averaging over initial conditions and time. To generate 
a non-zero dc current it is necessary to break all relevant symmetries which are respon- 
sible for transforming a trajectory with x(t; Xo,Po), p(t;xo,Po) , x(to;xo,po) = xo and 
p(to] xo,po) = po into another one with a momentum of the opposite sign. Note, that the 
generated trajectory must belong to the same statistical ensemble. These symmetries can 
be expressed via combinations of shifts and reflections in time and space (see 0): 



S a 


x(t; x ,p ) 
p(t;x ,p ) 




-x (t + T/2;x ,p ) + 2X 
-p (t + T/2;x ,p ) 


(2) 


s b 


x(t; x ,p ) 
p(t;x ,p ) 




x(-t + 2r; x ,p ) 
-p(-t + 2t;x ,Po) 




(3) 



Here the constants X and r are defined by the shape of functions U(x) and E(t) and are in 
correspondence with the discrete symmetry groups of the Hamiltonian ([I]). If U'(x + X) — 
—U\—x + X) and E(t) = —E(t + T/2) the dynamical equations of motion are invariant 
under the symmetry operation S a (this symmetry has been also observed and discussed by 
Ajdari et al for the overdamped case). On the other hand, if E(t + r) = E(—t + r), they 
have Sb symmetry (dissipationless case provided). Both symmetries (@-0) m ay be broken 
by choosing appropriate functions U(x) and E(t). 

In the overdamped case the corresponding dynamical equation 

7 x + U'(x) + E(t) = , (4) 

may possess an additional symmetry 
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x(t;x ,p ) 




pit] X ,Po) 





(5) 



x {—t + r; x ,p ) + X 
-p(-t + t;x ,Po) 

provided that U(x) and E(t) possess the following properties: U'(x) = —U'(x + X) and 
E(t) = —E(—t + r). Note that the symmetry S c if realized is based on time reversal in the 
overdamped limit, which is quite unexpected. An independent discovery of this symmetry 
has been done in (coined "supersymmetry" ) . 

The model analysis in [|| does not account for more realistic statistical properties of an 
ensemble of particles (like thermalization of the distribution). Still it provides one with a 
good intuitive understanding of current rectification. On the other hand, a lot of solid state 
applications require a more rigorous statistical description of transport properties. The aim 
of the present paper is to show that the symmetry approach can be applied as well to the 
classical kinetic Boltzmann equation. We will show in particular that the attractor of the 
Boltzmann equation completely reflects the symmetries of the original equation of motion 
for one particle. The kinetic Boltzmann equation reads (see, for example, ||13|| ) 



Cf = d t f + xd x f+pd p f = J(f,F) 
x = p , p — —U'(x) + E{t) . 



(6) 
(7) 



Here / = f(t,x,p) is an unknown distribution function, F = F(x,p) is some equilibrium 
distribution function chosen as: 



F(x,p) =F x {x)F p {p) 



e -p 2 /2 e -U(x) 



L= / e~ u{x) dx 
Jo 



(8) 
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and J{fi F) is a collision integral. Here we apply the momentum independent single re- 



laxation time approximation (also known as "r-approximation" ) [13|. Thus, the collision 
integral should be written as 



(9) 



where v is the dissipation constant or the characteristic relaxation frequency (1/V is the re- 
laxation time to the equilibrium state). This is the simplest form of a collision integral which 



includes relaxation to the equilibrium distribution and breaks time-reversal symmetry. We 
do not intend to discuss its validity for specific physical realizations. Rather we consider the 
resulting kinetic equation as a treatable case which allows us to demonstrate the qualitative 
features of dynamical symmetry breaking. Note that for any finite v time reversal symmetry 
is broken, so that we are concerned with S a only. The limit of v — > will lead to a restoring 
of time reversal symmetry (see below) and consequently to a necessary consideration of Sb 
as well. Interestingly the opposite limit v — > oo will also induce a symmetry S c . 

We start with a perturbative description of the overdamped case when \ jv is the smallest 
time scale in the system (formally, we put v ^> 1). We rewrite Eq. (||) in the following form 

f(x,p,t) = F(x,p)-±-Cf(x,p,t). (10) 

Since 1/u is a small parameter here, we expand the unknown function / into a series in 
terms of \jv: 

f(x, P ,t) = F(x, P ) + £(-iy^f^ . (ii) 

i=l v 

Collecting terms of the same order, we successively find the increments fi. The average 
value of the dc current is calculated using 

jdc = (j(t))t , j(t) = (pf(t,x,p)) x>p , (12) 

where {■■■) xpt stands for averaging over the phase space coordinates (x,p) and time cor- 
respondingly. We expect the attractor of the Boltzmann equation to possess the same 
symmetries as the corresponding dynamical equations of motion, which should be reflected 
in the vanishing or nonvanishing of jdc- In other words, if the relevant symmetries are broken 
in the equations of motion ||, we will observe this fact by analysing the increments f\. We 
show this using the following model 

U(x) = U jl - cos x + y [1 - cos(2x + 0)] j , (13) 
E{t) = Ex cos tut + E 2 cos(2wt + a) . (14) 
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As it was already shown in the previous work ||, this choice of the potential and external 
force can break both symmetries (|2])-(|3]). We apply the perturbation analysis for two specific 
cases: (i) v 2 7^ 0, E 2 = and (ii) v 2 = 0, E 2 ^ 0. A related overdamped case was discussed 



in terms of harmonic mixing in JT4|. In both cases the appropriate time-space symmetries 
are broken except for specific values of 9 and a. As already explained, we need to break 
only S a since time reversal and thus Sb are already violated for nonzero v. We find (a more 
detailed description of these calculations will be reported elsewhere) that the lowest order 
increment which yields a non-zero dc current in case (i) is proportional to 1/V 5 and to the 
square of the external field E(t): 

3 - r2n 



Jdc 



4 (E\t)) (U'(x)F^x)) x = ——-—rEf ru"\x)e- u ^dx . (15) 

V * 11 AU Lj JO 



The integral in the above formula cannot be computed analytically. Yet we find that it is 
non-zero only when v 2 ^ and 9 ^ 0, tt because otherwise the symmetry of the potential 
function enforces the vanishing of the above expression. For the case (ii) this term equals 
zero and the first non-zero contribution to the current appears in seventh order of the 
perturbation theory. As a result we obtain the following expression for the dc current 



J dc = ~ {E 3 (t)) t {U'{x)F-{x)) x = -^TJ^^cosa ■ (16) 



Here I n (z) is the modified Bessel function of nth. order [15|. Note that the averaged current 
is proportional to cos a in this case. The symmetry which causes the vanishing of the dc 
current for a = ±n/2 is S c [see Eq. (||)]. The above considered case (ii) with a = ±n/2 
satisfies exactly this symmetry. To conclude the consideration of ^ > 1 we find perfect 
agreement between the symmetry properties of the dynamical equations of motion and the 
symmetry properties of the used kinetic equation, resulting in a vanishing or appearance of 
a dc current. 

Let us now investigate the case of intermediate v ~ 1 and weak v <^ 1 dissipation. 
When the characteristic relaxation frequency v approaches zero the solution of the kinetic 
equation is expected to adequately describe the properties of the (microcanonical) ensemble 



of trajectories in the phase space of the underlying dynamical system. Especially we expect 
that if allowed by the choice of functions U(x) and E(t), the symmetry Sb described in 
detail in || should be restored on the attractor of the kinetic equation in the limit of v — > 0. 

The limit of small dissipation appears to be difficult to be treated analytically. We 
have used numerical simulations to calculate the dc current in this case. Two independent 
numerical methods have been used. 

The first one (method I) is based on the expansion of the unknown distribution function 
f(t,x,p) into Hermitian polynomials in p and Fourier series in x-space. First we expand / 
into a series w.r.t. p: 

f{t,z,p) = j£C n {t,z)\n) , |n)= J^-. e-^ (17) 
n =o 7r i / 4 v2 n n! 

where H n (p) are Hermite polynomials of n-th order. Since U(x) and F(x,p) are periodic 
functions in x we may use 

\r\\ 00 00 00 

F ^p) = wrrr± £ F ^ imx > u ( x ) = E u m e im * , c n (t,x)= Km(ty mx . 

^ m=—oo m=—oo m=—oo 

(18) 

Using (n\n') = 5 n y we arrive at the following infinite set of coupled ordinary differential 
equations for n = 

] ^ 1 

A , m + -^=E(t)A ljm - ^ /f/zA n>m _ ; + -^=my4i jm = -v{A , m - F m ) (19) 

and n > 



Ai,m + -j=E(t)(Vn + L4 n +i,m - v^^n-i.m) + — 7=m(Vn + 1 A n+ i jOT + VnAi-i,r 



V2 1 \J I { 

The resulting ODEs for the time-dependent coefficients A njm (t) are integrated numerically. 
After a sufficiently large integration time ^> 1/v the system will relax into a time-periodic 
attractor. On this attractor the induced current is given by 
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oo 



f(2nY 

( p ) P:X = 3^2^ ^2„ + i,o^ r y : V2^Tl (21) 

n=0 ^ U - 

which can be further averaged over time. Actual computation requires cutoffs in n and m. 
Results presented below are checked to be independent on the chosen cutoff values. For the 
following results the cutoff values were n max = 60 and m max = 30. 

The second approach (method II) is based on the method of characteristics. Taken an 
equilibrium distribution at zero time, the solution of the linear kinetic equation with the 
constant relaxation time can be expressed through the equilibrium distribution function and 



a statistically weighted contribution of the trajectories in the phase space ||16|| : 
f(x, P ,t)=F (t, t , x, p), r)!f\t, t , x, p) 



e -K*-*o) _|_ 



+ v re^ {t - t ' ) Fk(t / ,x,p),r/ p (t',x,p)ldt / . (22) 

J in 



'to 

We remind that the functions r]^ stand for the inverted equation of a trajectory, i.e., 
they denote the dependence of the point on a trajectory {x,p, t} on its initial conditions; 
characteristics r] X)P correspond to the trajectory which passes through the point {x,p, t'} 
in the phase space. For brevity we put also to = 0. Since the underlying dynamical system 
is non-integrable we have to find the characteristics numerically. We choose a large random 
ensemble of initial conditions with \p Q \ < \p ma x\ and use them directly to generate an 
approximation of f(x,p,t) in accordance with Eq.(^). The accuracy of the approximation 
is controlled by the number of trajectories (initial conditions) used. We used up to 10 5 
trajectories to gain good averaging statistics. The net error of calculations was estimated as 
~ 10 — 20% for large relaxation parameter v or large dc current jd c after averaging of the 
current j(t) over several periods T of the field E(t). The maximal error has been detected in 
those cases when the average dc current tends to zero, for example, when we put v = 0.001 
and E 2 = (see inset of Fig. 2). In this case the fluctuations of the mean value of the 
current averaged over several periods T were of the order of residual value of the dc current 
(Aj dc ~ \j dc \ = 0.002 at v = 0.001 and E 2 = 0). In order to achieve the accuracy of the order 
of Ajdc/jdc ~ 1/3 -i- 1/2 we had to perform averaging over hundreds of periods T. Note that 



method II has some analogy to the direct sampling in classical Monte-Carlo simulations. 

Method I is an optimal tool for not too small relaxation v > 0.01 while method II shows 
up with larger errors in calculations but covers, nevertheless, the limit of v —>■ 10~ 3 . Both 
methods give nearly identical results for v > 10~ 2 -v- 10 _1 . 

Fig. [I] displays the dependence of the full time-dependent current on time on the attractor 
of the kinetic equation computed with the help of method II for v = 0.01. We observe the 
expected periodic variation in time, with a nonzero time averaged dc current indicated by 
the thick horizontal line. Note that in this case the dc current is of the order of the time- 
dependent current amplitude (only three times less). This implies that the underdamped 
case may serve as an ideal testing ground in various applications, since the rectification effect 
is of the order of 30 percent. 
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FIG. 1. Time dependence of j(t) on the attractor of the kinetic equation (thin solid line). The 
thick solid horizontal line shows the value of the dc part of j(t). Parameters: v = 0.01, a = 1.5, 
V2 = 0, Uq = 6, Ei = —2.6, E<i = —2.04, u = 0.85. The dashed line denotes the reference of j = 0. 

In Fig. we show the dependence of the averaged dc current on the dissipation parameter 
v. Results from both above described methods are presented. They show data which are 
similar qualitatively and quantitatively (see the caption for details). 
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FIG. 2. Dependence of the averaged dc current (|T^), on the dissipation parameter u, computed 
for the case of broken time symmetry using numerical methods I (o) and II (*) with parameters 
as in Fig.|l] except a = 1. The inset shows a log-log plot of this dependence with parameters as in 
Figjl| except a = (o) and E 2 = 0, v 2 = 0.8, 9 = 2 (+). 

Consider the case of E 2 ^ 0, a ^ 0, v 2 = 0. Directed current appears due to a violation 
of the S a symmetry. When the dissipation v reaches zero, this symmetry is still broken, 
however, the dynamical equations of motion may restore a second symmetry, Sb- But for 
the above set of parameters, Sb is also violated, and a non-zero dc current should persist 
down to zero dissipation. This is observed in Fig. 0. Suppose now we choose a = 0, then S a 
is still violated, and jd c ^ for non-zero dissipation (see inset of Fig. g). However, in the 
limit of zero dissipation Sb is restored: E(t) = E(—t) and our numerical simulations show 
that current tends to zero as predicted by symmetry considerations. The same happens in the 
case of a "ratchet potential", v 2 ^ and without harmonic mixing (E 2 = 0). For v ^ both 
S a and Sb are violated, resulting in a nonzero dc current. Again the limit of zero dissipation 
leads to a restoring of Sb- This leads to a disappearance of the average dc current at v = 
(see inset of Fig. 0). Thus the numerical simulations confirm our expectation that the 
attractor properties of the kinetic equation reflect the symmetries S a , Sb of the underlying 
dynamical equations for all v including v — > 0. 
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In Fig. [3] we show the dependence of the averaged dc current on the "phase difference" 
a for three different values of the dissipation parameter v = 4,1,0.3, computed with the 
help of method I. Note that the currents for v — 4 and v = 1 are scaled with factors 181.6 
and 4.86 respectively, consequently we observe an enhancement of the maximum dc current 
value by more than two orders of magnitude when going from the overdamped case [y = 4) 
to the underdamped one [y = 0.3). Besides this observation Fig. fj displays again the 
consequences of symmetry breaking or restoral. With increasing v the dc current becomes 
zero for a = ±7r/2 which is a consequence of restoring of S c symmetry. With decreasing v 
we observe a vanishing of the current at a = 0, tt which is a consequence of restoring of the 
Sb symmetry. 

1 i ■ , ■ , , , , 1 
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FIG. 3. Dependence of the averaged dc current on the phase delay parameter a for different 
values of v: u = 0.3 (filled circles), v = 1 (squares) and v = 4 (triangles). Note that the current 
values are scaled by a factor of 4.86 (y = 1) and 181.6 (y = 4). The lines connecting the data are 
splines. Other parameters are as in Fig.|l]. 

Conclusions. In this paper, we extend the symmetry approach for explanation of a dc 
current in driven nonlinear systems to the framework of the classical kinetic Boltzmann 
equation. We have analysed separately the cases of large and small dissipation and have 
demonstrated that the rectification of current due to nonlinear mixing of harmonics takes 
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place in agreement with previous symmetry consideration ||. In particular, it follows that 
in the limit of small v the attractor of the Boltzmann equation (|6|-[7|) reflects the dynamical 
symmetries of the Newtonian equations of motion [[T7|]. Finally we observe that a proper 
choice of the system parameters in the underdamped case ensures the value of dc current to 
be of the order of the ac current amplitude. The underdamped case leads to an enhancement 
of the dc current by several orders of magnitude as compared to the overdamped one. 
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and stimulating discussions. We also would like to thank P. Reimann for sending us preprints 
prior to publication. 
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